Ahora queremos estudiar cómo encontrar los parámetros óptimos de un modelo para minimizar una función de pérdida. Existen muchos algoritmos de optimización. Uno de los que converge en menos pasos es el método de Newton, pero su costo computacional crece rápidamente con el número de parámetros. Las redes neuronales pueden involucrar billones de parámetros, lo que hace que el método de Newton y otros similares sean imposibles de aplicar. Por lo tanto se usa el Descenso de Gradiente (Gradient Descent, GD) al entrenar redes neuronales.
El descenso de gradiente consiste en tomar el gradiente de la función de pérdida en el espacio de parámetros. La idea es ir en la dirección opuesta al gradiente, ya que el gradiente apunta en la dirección de máximo crecimiento de la función (y su opuesto en la dirección de máximo descenso). Cuando llegamos a un mínimo de la función, el gradiente será cero.
En esta sección veremos primero cómo aplicar el descenso de gradiente para un perceptrón multicapa asumiendo que el gradiente de la función de pérdida con respecto a los parámetros de la red es conocido. Luego veremos cómo calcular el gradiente de la función de pérdida con respecto a los parámetros de la red de manera eficiente. Finalmente discutiremos algunos detalles de cómo escoger el punto inicial en el espacio de parámetros para inciar el descenso de gradiente.
Descenso de Gradiente
Ejemplo de Regresión Lineal Simple con Descenso de Gradiente
Antes de aplicar el descenso de gradiente a una red neuronal, empezamos por hacerlo en un caso sencillo que es fácil de visualizar: la regresión lineal simple. En este caso queremos ajustar la pendiente y el intercepto de una recta \(y = m x + b\) a un conjunto de datos \((x_i, y_i)\). La función de pérdida (MSE) a minimizar es: \[
\mathcal{J}(m,b) = \frac{1}{N}\sum_{i=1}^N (y_i - (m x_i + b))^2.
\] donde \(N\) es el número total de puntos de datos. Notemos que la función de pérdida es una función de dos variables, \(m\) y \(b\). Por eso la podemos graficar fácilmente.
Código
import numpy as npimport matplotlib.pyplot as pltfrom mpl_toolkits.mplot3d import Axes3D# Generar muestra aleatorianp.random.seed(42)N =100x = np.random.uniform(-10, 10, N)m_verdadero, b_verdadero =2, -1y = m_verdadero * x + b_verdadero + np.random.normal(0, 2, N)# Definir función de pérdidadef perdida_ecm(m, b, x, y):return np.mean((y - (m * x + b))**2)# Crear malla para gráfico 3D y heatmaprango_m = np.linspace(-3, 7, 100)rango_b = np.linspace(-24, 24, 100)M, B = np.meshgrid(rango_m, rango_b)Z = np.array([perdida_ecm(m, b, x, y) for m, b inzip(M.ravel(), B.ravel())]).reshape(M.shape)# Configurar el layout de los subplotsfig = plt.figure(figsize=(6, 12))# Graficar superficie 3Dax1 = fig.add_subplot(211, projection='3d')superficie = ax1.plot_surface(M, B, Z, cmap='viridis')ax1.set_xlabel('m (pendiente)')ax1.set_ylabel('b (intercepto)')ax1.set_zlabel('Pérdida ECM')ax1.set_title('Función de Pérdida 3D')# Graficar heatmap 2Dax2 = fig.add_subplot(212)heatmap = ax2.imshow(Z, extent=[rango_m.min(), rango_m.max(), rango_b.min(), rango_b.max()], origin='lower', cmap='viridis', aspect='auto')contours = ax2.contour(M, B, Z, colors='black', levels=10)ax2.set_xlabel('m (pendiente)')ax2.set_ylabel('b (intercepto)')ax2.set_title('Heatmap de la Función de Pérdida')plt.tight_layout()plt.show()
Figura 1: Función de pérdida para regresión lineal
En la Figura 1 se muestra la función de pérdida para regresión lineal. Vemos que esta es mínima en el punto que corresponde a los valores verdaderos de \(m\) y \(b\) y que es convexa, es decir, que cualquier línea recta que conecte dos puntos está por encima de la superficie.
Utilizamos el Descenso de Gradiente (Gradient Descent, GD) para encontrar los valores óptimos de \(m\) y \(b\). La actualización en cada iteración se realiza siguiendo la dirección opuesta al gradiente de la función de pérdida:
\[
\begin{aligned}
m &\leftarrow m - \eta \,\frac{\partial \mathcal{L}}{\partial m},\\
b &\leftarrow b - \eta \,\frac{\partial \mathcal{L}}{\partial b},
\end{aligned}
\]
Aquí, \(\eta\) es la tasa de aprendizaje (learning rate), un hiperparámetro que controla el tamaño de los pasos que damos en cada actualización. Las derivadas parciales (el gradiente) indican la dirección de máximo ascenso de la pérdida, y su inverso es la dirección de máximo descenso. Al restar el gradiente multiplicado por \(\eta\), nos movemos un poco hacia el mínimo. En este caso podemos calcular fácilmente el gradiente de la función de pérdida con respecto a \(m\) y \(b\): \[
\frac{\partial \mathcal{J}}{\partial m}
= -\frac{2}{N}\sum_{i=1}^N (y_i - (m x_i + b))\,x_i,
\]\[
\frac{\partial \mathcal{J}}{\partial b}
= -\frac{2}{N}\sum_{i=1}^N (y_i - (m x_i + b)).
\]
Código
import matplotlib.animation as animationfrom IPython.display import HTML# Definir función de pérdida y sus derivadasdef perdida_ecm(m, b, x, y):return np.mean((y - (m * x + b))**2)def derivada_m(m, b, x, y):return-2* np.mean((y - (m * x + b)) * x)def derivada_b(m, b, x, y):return-2* np.mean(y - (m * x + b))# Crear malla para gráfico 3D y heatmaprango_m = np.linspace(0, 4, 100)rango_b = np.linspace(-5.5, 2, 100)M, B = np.meshgrid(rango_m, rango_b)Z = np.array([perdida_ecm(m, b, x, y) for m, b inzip(M.ravel(), B.ravel())]).reshape(M.shape)# Configurar figurafig, ax = plt.subplots(figsize=(10, 6))heatmap = ax.imshow(Z, extent=[rango_m.min(), rango_m.max(), rango_b.min(), rango_b.max()], origin='lower', cmap='viridis', aspect='auto')contours = ax.contour(M, B, Z, colors='black', levels=10)ax.set_xlabel('m (pendiente)')ax.set_ylabel('b (intercepto)')ax.set_title('Descenso de Gradiente')# Ejecutar el descenso de gradientem_inicio, b_inicio =3, -2.5# Punto inicialeta =0.015# Tasa de aprendizajepasos =15# Número de pasos a mostrar# Calcular la trayectoria completa de antemanotrayectoria_m = [m_inicio]trayectoria_b = [b_inicio]m_actual, b_actual = m_inicio, b_iniciofor _ inrange(pasos): grad_m = derivada_m(m_actual, b_actual, x, y) grad_b = derivada_b(m_actual, b_actual, x, y)# Actualizar parámetros m_actual -= eta * grad_m b_actual -= eta * grad_b# Guardar trayectoria trayectoria_m.append(m_actual) trayectoria_b.append(b_actual)# Inicializar elementos de la animaciónline, = ax.plot([], [], 'ro-', linewidth=2, markersize=8)points = ax.scatter([], [], color='red', s=50)# Función de inicializacióndef init(): line.set_data([], []) points.set_offsets(np.empty((0, 2)))return line, points# Función de actualización para cada framedef update(frame): line.set_data(trayectoria_m[:frame+1], trayectoria_b[:frame+1]) points.set_offsets(np.column_stack((trayectoria_m[:frame+1], trayectoria_b[:frame+1])))return line, points# Crear animaciónani = animation.FuncAnimation(fig, update, frames=len(trayectoria_m), init_func=init, blit=True, interval=300)HTML(ani.to_jshtml())
(a) Descenso de gradiente para regresión lineal
(b)
Figura 2
En la Figura 2 se muestra el descenso de gradiente para regresión lineal. Graficamos unos pocos pasos del descenso de gradiente. Ya podemos apreciar un posible problema: El descenso avanza rápidamente en la dirección de la pendiente, pero luego se detiene y avanza lentamente en la dirección del intercepto. Esto es porque el gradiente es más pronunciado en la dirección de la pendiente que en la dirección del intercepto. Veremos más adelante cómo intentar mejorar esto.
Ejemplo 2D de modelo Gabor: paisajes con mínimos locales
El ejemplo de arriba fue para una función de pérdida convexa. En general las funciones de pérdida de las redes neuronales no tienen esta propiedad. Por lo tanto, el “paisaje” de la pérdida puede ser complejo, presentando múltiples mínimos locales y puntos de silla (saddle points). Esto dificulta la optimización, ya que el Descenso de Gradiente puede quedar atrapado en un mínimo local que no es el mínimo global. Usaremos una función sencilla con esas características llamada función de Gabor para ilustrar el problema.
Usamos una función Gabor 1D como modelo, donde ajustaremos dos de sus parámetros. Consideramos \(\sigma\) y \(\omega\) como parámetros fijos conocidos. Optimizaremos la amplitud \(A\) y la fase \(\phi\): \[
f(x; \phi_0,\phi_1) = e^{-\frac{(\phi_0 + \phi_1 x)^2}{32}}\cos(\phi_0 + \phi_1 x + \phi).
\] Ajustamos \((\phi_0,\phi_1)\) a datos \((x_i, y_i)\) generados usando valores verdaderos \((\phi_0^*,\phi_1^*)\), en este ejemplo.
Minimizamos el Error Cuadrático Medio entre las predicciones del modelo Gabor y los datos \(y_i\). Siendo \(N\) el número de puntos \(x_i\): \[
\mathcal{L}(\phi_0,\phi_1) = \frac{1}{N}\sum_{i=1}^N (y_i - f(x_i; \phi_0,\phi_1))^2.
\] La naturaleza oscilatoria del término \(\cos(\phi_0 + \phi_1 x + \phi)\), multiplicada por la amplitud \(A\), crea un paisaje de pérdida no convexo con múltiples valles (mínimos locales) y puntos de silla.
Para entender el paisaje de la pérdida, generamos un gráfico de contorno de \(\mathcal{L}(\phi_0,\phi_1)\) evaluando la pérdida en una rejilla de valores de \(\phi_0\) y \(\phi_1\). Este gráfico muestra las regiones de baja pérdida (valles) y las zonas de puntos de silla. Adicionalmente, se pueden trazar las trayectorias seguidas por el algoritmo de Descenso de Gradiente partiendo desde diferentes puntos iniciales \((\phi_{0,0}, \phi_{1,0})\).
Código
# Definir función de Gabor y generar datosdef gabor(x, phi_0, phi_1, sigma=5.0):return np.exp(-(phi_0 + phi_1*x)**2/(2*sigma**2)) * np.sin(phi_0 + phi_1*x)def d_gabor_dx(x, sigma=5.0):return- x * np.exp(-x**2/(2*sigma**2))/sigma**2* np.sin(x) +\ np.exp(-x**2/(2*sigma**2)) * np.cos(x)# Generar datos sintéticosnp.random.seed(42)N =100x = np.linspace(-3, 3, N)phi_0_verdadero, phi_1_verdadero =2.0, 1.0sigma =5.0# Parámetros fijosy = gabor(x, phi_0_verdadero, phi_1_verdadero, sigma) + np.random.normal(0, 1, N)# Definir función de pérdida y sus derivadasdef perdida_ecm(phi_0, phi_1, x, y, sigma=5.0): y_pred = gabor(x, phi_0, phi_1, sigma)return np.mean((y - y_pred)**2)def derivada_phi_0(phi_0, phi_1, x, y, sigma=5.0): y_pred = gabor(x, phi_0, phi_1, sigma)return-2*np.mean((y - y_pred) * d_gabor_dx(phi_0 + phi_1*x, sigma))def derivada_phi_1(phi_0, phi_1, x, y, sigma=5.0): y_pred = gabor(x, phi_0, phi_1, sigma)return-2*np.mean((y - y_pred) * d_gabor_dx(phi_0 + phi_1*x, sigma)*x)# Crear malla para gráfico de contornorango_phi_0 = np.linspace(-1, 4, 100)rango_phi_1 = np.linspace(-2, 4, 100)Phi_0_grid, Phi_1_grid = np.meshgrid(rango_phi_0, rango_phi_1)Z = np.array([perdida_ecm(p_0, p_1, x, y) for p_0, p_1 inzip(Phi_0_grid.ravel(), Phi_1_grid.ravel())]).reshape(Phi_0_grid.shape)# Configurar figurafig, ax = plt.subplots(figsize=(10, 6))heatmap = ax.imshow(Z, extent=[rango_phi_0.min(), rango_phi_0.max(), rango_phi_1.min(), rango_phi_1.max()], origin='lower', cmap='viridis', aspect='auto')contours = ax.contour(Phi_0_grid, Phi_1_grid, Z, colors='black', levels=10)ax.set_xlabel(r'$\phi_0$')ax.set_ylabel(r'$\phi_1$')ax.set_title('Descenso de Gradiente en Modelo Gabor')# Definir puntos iniciales para diferentes trayectoriaspuntos_iniciales = [ (2.3, 1.0), # (φ_0, φ_1) - Punto 1 (1.0, 3.0), # Punto 2 (3.0, 0.0), # Punto 3 (1.0, -1.0), # Punto 4 (2.5, -1.5) # Punto 5]colores = ['r', 'b', 'g', 'y', 'c']eta =0.1# Tasa de aprendizajepasos =30# Número de pasos a mostrar# Calcular todas las trayectorias de antemanotodas_trayectorias = []for phi_0_inicio, phi_1_inicio in puntos_iniciales: trayectoria_phi_0 = [phi_0_inicio] trayectoria_phi_1 = [phi_1_inicio] phi_0_actual, phi_1_actual = phi_0_inicio, phi_1_iniciofor _ inrange(pasos): grad_phi_0 = derivada_phi_0(phi_0_actual, phi_1_actual, x, y) grad_phi_1 = derivada_phi_1(phi_0_actual, phi_1_actual, x, y)# Actualizar parámetros phi_0_actual -= eta * grad_phi_0 phi_1_actual -= eta * grad_phi_1# Guardar trayectoria trayectoria_phi_0.append(phi_0_actual) trayectoria_phi_1.append(phi_1_actual) todas_trayectorias.append((trayectoria_phi_0, trayectoria_phi_1))# Inicializar elementos de la animaciónlines = []points_collection = []for i, color inenumerate(colores): line, = ax.plot([], [], f'{color}-', linewidth=2, markersize=8) points = ax.scatter([], [], color=color, s=50) lines.append(line) points_collection.append(points)# Función de inicializacióndef init():for line in lines: line.set_data([], [])for points in points_collection: points.set_offsets(np.empty((0, 2)))return lines + points_collection# Función de actualización para cada framedef update(frame):for i, ((trayectoria_phi_0, trayectoria_phi_1), line, points) inenumerate(zip(todas_trayectorias, lines, points_collection)): line.set_data(trayectoria_phi_0[:frame+1], trayectoria_phi_1[:frame+1]) points.set_offsets(np.column_stack((trayectoria_phi_0[:frame+1], trayectoria_phi_1[:frame+1])))return lines + points_collection# Crear animaciónani = animation.FuncAnimation(fig, update, frames=pasos+1, init_func=init, blit=True, interval=300)HTML(ani.to_jshtml())
(a) Descenso de gradiente para modelo Gabor
(b)
Figura 3
En la Figura 3 se muestra el paisaje de la pérdida para el modelo Gabor, y las trayectorias de convergencia para diferentes puntos iniciales. Vemos que diferentes puntos van hacia diferentes mínimos locales, lo que dificulta la convergencia al mínimo global.
Descenso de Gradiente Estocástico (SGD) y por mini-lotes
Cuando el conjunto de datos es muy grande, calcular el gradiente sobre todos los \(N\) ejemplos en cada paso puede ser computacionalmente muy costoso. Por ello, se utilizan variantes del GD que usan un subconjunto de los datos para aproximar el gradiente a cada paso. Esta será sólo una aproximación, y nos da una estima ruidosa del gradiente de la pérdida ideal. Pero esto tiene una gran ventaja: Ese ruido puede empujar el algoritmo a saltar de mínimos locales poco profundos. Veamos dos variaciones de esa idea:
Descenso de Gradiente Estocástico (SGD) Realiza una actualización de parámetros por cada ejemplo individual \((x_i, y_i)\). En cada paso, el gradiente se estima usando solo ese único ejemplo: \(\theta \leftarrow \theta - \eta \nabla_{\theta} \mathcal{L}_{\text{indiv}}(\theta; x_i, y_i)\), donde \(\mathcal{L}_{\text{indiv}}\) es la pérdida calculada para ese punto. Luego repite para todos los puntos \((x_i, y_i)\) y vuelve a empezar.
Ventajas: Actualizaciones muy rápidas y bajo coste computacional por actualización. La alta varianza en la estimación del gradiente puede ayudar a escapar de mínimos locales “malos”.
Desventajas: Al usar un único dato para calcular el gradiente, la trayectoria de convergencia es muy ruidosa (alta varianza). Puede que nunca converja exactamente al mínimo, sino que oscile alrededor de este.
Descenso de Gradiente por Mini-lotes: Es un compromiso entre los dos anteriores y el método más común en Deep Learning. El gradiente se calcula sobre un pequeño subconjunto aleatorio de datos llamado mini-batch (lote), de tamaño \(B\) (e.g., \(B=32, 64, 128\)). Luego repite para todos los lotes y vuelve a empezar.
Ventajas: Reduce significativamente la varianza del gradiente comparado con SGD, llevando a una convergencia más estable. Mantiene una alta eficiencia computacional aprovechando la paralelización de hardware (GPUs).
Desventajas: Introduce un nuevo hiperparámetro (tamaño del batch, \(B\)).
En la práctica, el descenso de gradiente por mini-lotes es el método más común y eficiente. Es el usado para entrenar las redes neuronales famosas que han impulsado el interés en el aprendizaje automático.
A continuación el pseudocódigo para el descenso de gradiente por mini-lotes:
# Inicializar parámetrostheta = theta_0# Bucle de entrenamientofor epoch inrange(num_epochs):# Mezclar datos np.random.shuffle(data)# Bucle sobre lotesfor batch in data:# Calcular gradiente grad = compute_gradient(batch)# Actualizar parámetros theta = theta - eta * grad
Momentum
Otra estrategia para atacar el problema de los mínimos locales y la oscilación alrededor de mínimos, es cambiar el proceso de aprendizaje. Es decir, modificamos la forma como se actualizan los pesos a cada paso.
El método de Momentum introduce un término de “velocidad” \(v\) que acumula una media móvil exponencial de los gradientes pasados. Ayuda a acelerar el descenso en direcciones consistentes y amortigua las oscilaciones en direcciones donde el gradiente cambia rápidamente (típico en valles estrechos y alargados).
La actualización se modifica como sigue (donde \(\theta\) representa los parámetros del modelo, como \(m\) y \(b\)):
\(v_t\) es el vector de velocidad en el paso \(t\).
\(\beta\) es el coeficiente de momentum (usualmente cercano a 0.9). Controla cuánto del momentum anterior se conserva. Un \(\beta\) alto significa que los gradientes pasados tienen más influencia.
\(\nabla\mathcal{L}(\theta_{t-1})\) es el gradiente calculado en el paso actual (puede ser batch, mini-batch o stochastic).
A cantidades que se acumulan de la manera como lo hace \(v_t\) se les llama media móvil exponencial.
Esto es burdamente similar a una bola pesada rodando por una pendiente. Acumula momento, lo que le permite pasar por pequeños baches y acelerar en las bajadas consistentes. A pesar de llamarse momentum, no es exactamente el momentum definido en la física.
A continuación el pseudocódigo para el descenso de gradiente con momentum en el problema de regresión lineal:
# Inicializar velocidadesv_m, v_b =0.0, 0.0beta =0.9# Coeficiente de momentumeta =0.01# Tasa de aprendizaje# Bucle de entrenamiento (puede ser sobre épocas y mini-batches)for iteration inrange(num_iterations):# Calcular gradientes (usando batch, mini-batch, o SGD)# grad_m = ∂L/∂m, grad_b = ∂L/∂b grad_m, grad_b = compute_gradient(m, b, data_batch)# Actualizar velocidades v_m = beta * v_m + (1- beta) * grad_m v_b = beta * v_b + (1- beta) * grad_b# Actualizar parámetros usando las velocidades m -= eta * v_m b -= eta * v_b
Adam (Adaptive Moment Estimation)
Cuando hicimos el ejemplo de regresión lineal, vimos que la convergencia hacia el mínimo en una dimensión era mucho más rápida que en la otra. Adam es un optimizador que intenta resolver este problema ajustando automáticamente la tasa de aprendizaje para cada parámetro individualmente.
Adam es un optimizador muy popular que combina dos ideas principales:
Momentum: Utiliza una media móvil exponencial del gradiente (primer momento, \(m_t\)), similar al método de Momentum.
Adaptación de Tasa de Aprendizaje: Utiliza una media móvil exponencial del cuadrado del gradiente (segundo momento, \(v_t\)). Esto reduce la diferencia entre las tasas de aprendizaje más rápidas y las más lentas.
Además, aplica una corrección del sesgo (bias correction) a ambas estimaciones de momentos para contrarrestar su tendencia a estar sesgadas hacia cero durante las primeras iteraciones.
Las ecuaciones de actualización son: \[
\begin{aligned}
m_t &= \beta_1 m_{t-1} + (1-\beta_1)\,g_t & \text{(Estimación del 1er momento)} \\
v_t &= \beta_2 v_{t-1} + (1-\beta_2)\,g_t^2 & \text{(Estimación del 2do momento)} \\
\hat m_t &= \frac{m_t}{1-\beta_1^t} & \text{(Corrección del sesgo 1er momento)} \\
\hat v_t &= \frac{v_t}{1-\beta_2^t} & \text{(Corrección del sesgo 2do momento)} \\
\theta_t &\leftarrow \theta_{t-1} - \eta\,\frac{\hat m_t}{\sqrt{\hat v_t}+\epsilon} & \text{(Actualización del parámetro } \theta \text{)}
\end{aligned}
\]
Donde: * \(t\) es el número de iteración (paso de actualización). * \(g_t = \nabla_{\theta} \mathcal{L}(\theta_{t-1})\) es el gradiente en el paso \(t\). * \(\beta_1\) es el coeficiente de decaimiento para el primer momento (usualmente ~0.9). * \(\beta_2\) es el coeficiente de decaimiento para el segundo momento (usualmente ~0.999). * \(\eta\) es la tasa de aprendizaje global (e.g., 1e-3). * \(\epsilon\) es una constante pequeña (usualmente ~1e-8) para evitar la división por cero y mejorar la estabilidad numérica.
A continuación se muestra un pseudocódigo para el descenso de gradiente con Adam en el problema de regresión lineal.
# Inicializar momentosm_m, m_b =0.0, 0.0# Primer momento (media móvil del gradiente)v_m, v_b =0.0, 0.0# Segundo momento (media móvil del cuadrado del gradiente)beta1 =0.9# Coeficiente de decaimiento para el primer momentobeta2 =0.999# Coeficiente de decaimiento para el segundo momentoepsilon =1e-8# Constante para evitar división por ceroeta =0.01# Tasa de aprendizajet =0# Contador de iteraciones# Bucle de entrenamientofor iteration inrange(num_iterations): t +=1# Calcular gradientes (usando batch, mini-batch, o SGD) grad_m, grad_b = compute_gradient(m, b, data_batch)# Actualizar primer momento (media móvil del gradiente) m_m = beta1 * m_m + (1- beta1) * grad_m m_b = beta1 * m_b + (1- beta1) * grad_b# Actualizar segundo momento (media móvil del cuadrado del gradiente) v_m = beta2 * v_m + (1- beta2) * (grad_m**2) v_b = beta2 * v_b + (1- beta2) * (grad_b**2)# Corregir el sesgo m_m_corrected = m_m / (1- beta1**t) m_b_corrected = m_b / (1- beta1**t) v_m_corrected = v_m / (1- beta2**t) v_b_corrected = v_b / (1- beta2**t)# Actualizar parámetros m -= eta * m_m_corrected / (np.sqrt(v_m_corrected) + epsilon) b -= eta * m_b_corrected / (np.sqrt(v_b_corrected) + epsilon)
Diferenciación Automática (AD)
Vimos que el descenso de gradiente y los algoritmos inspirados en este usan la derivada de la función de pérdida con respecto a los parámetros de la red. Las redes neuronales tienen hasta billones de parámetros, por lo que calcular la derivada de la función de pérdida con respecto a cada uno de ellos puede ser bastante costoso. Sin embargo, se han desarrollado técnicas que atacan este problema de forma bastante ingeniosa.
La diferenciación automática (AD) es un conjunto de técnicas que permite calcular de forma exacta y eficiente derivadas de funciones definidas por programas informáticos. Es crucial distinguirla de la diferenciación simbólica (que puede llevar a expresiones muy complejas) y de la diferenciación numérica (que usa aproximaciones por diferencias finitas y sufre de errores de truncamiento y redondeo). AD calcula derivadas exactas (hasta la precisión de la máquina) con una eficiencia comparable a la evaluación de la propia función.
Algoritmo de Backpropagation (Modo Reverso de AD)
En redes neuronales y otros modelos con muchos parámetros, la diferenciación automática se implementa de forma eficiente mediante el modo reverso (reverse mode), cuyo algoritmo más conocido es backpropagation. Existe también el modo directo (forward mode), que es menos eficiente pero puede ser útil en ciertas situaciones, y no lo veremos aquí por falta de tiempo.
El modo reverso calcula todas las derivadas de una única salida (escalar, como la función de pérdida \(\mathcal{J}\)) con respecto a todas las entradas y parámetros intermedios. Lo hace propagando las derivadas hacia atrás, desde la salida hacia la entrada, aplicando eficientemente la regla de la cadena sobre el grafo computacional definido por la evaluación de la función (el forward pass).
El proceso tiene dos fases:
Forward pass (Paso hacia adelante)
Partimos de la entrada \(x\) y calculamos, capa a capa, las activaciones intermedias y la salida final (y la pérdida \(\mathcal{J}\)). \[
z^{(1)} = W^{(1)} x + b^{(1)}, \quad a^{(1)} = \sigma(z^{(1)}), \quad z^{(2)} = W^{(2)} a^{(1)} + b^{(2)}, \dots, \quad \mathcal{J} = \text{Loss}(\text{output}, \text{target})
\]
Se almacenan en memoria todas las variables intermedias (\(z^{(\ell)}\), \(a^{(\ell)}\), \(W^{(\ell)}\), \(x\)) que influyen en el resultado final, ya que serán necesarias para calcular las derivadas en el paso inverso.
Backward pass (Paso hacia atrás)
Se inicializa con el gradiente de la pérdida respecto a sí misma (\(\frac{\partial \mathcal{J}}{\partial \mathcal{J}} = 1\)) o respecto a la salida final de la red antes de la función de pérdida.
Se propaga el gradiente hacia atrás capa por capa (\(\ell = L, L-1, \dots, 1\)). Para cada capa \(\ell\):
Se parte del gradiente respecto a la activación de la capa, \(\frac{\partial \mathcal{J}}{\partial a^{(\ell)}}\) (obtenido de la capa \(\ell+1\), o inicial si \(\ell=L\)).
Se calcula el gradiente respecto a la entrada lineal de la capa \(\ell\), denotado \(\delta^{(\ell)}\): \[
\delta^{(\ell)} = \frac{\partial \mathcal{J}}{\partial z^{(\ell)}} = \frac{\partial \mathcal{J}}{\partial a^{(\ell)}} \odot \sigma'(z^{(\ell)})
\] donde \(\odot\) es el producto elemento a elemento (Hadamard) y \(\sigma'\) es la derivada de la función de activación \(\sigma\) aplicada elemento a elemento a \(z^{(\ell)}\). En el caso de ReLU, \(\sigma'(z^{(\ell)}) = 1\) para \(z^{(\ell)} > 0\) y \(\sigma'(z^{(\ell)}) = 0\) para \(z^{(\ell)} \leq 0\).
Se calculan los gradientes respecto a los parámetros de la capa \(\ell\) usando \(\delta^{(\ell)}\) y la activación de la capa anterior \(a^{(\ell-1)}\) (guardada durante el forward pass): \[
\frac{\partial \mathcal{J}}{\partial W^{(\ell)}} = \delta^{(\ell)} \, (a^{(\ell-1)})^T
\]\[
\frac{\partial \mathcal{J}}{\partial b^{(\ell)}} = \delta^{(\ell)} \quad \text{(sumado sobre las muestras del batch si aplica)}
\]
Se calcula el gradiente que se propagará a la capa anterior \(\ell-1\): \[
\frac{\partial \mathcal{J}}{\partial a^{(\ell-1)}} = (W^{(\ell)})^T \, \delta^{(\ell)}
\] Este \(\frac{\partial \mathcal{J}}{\partial a^{(\ell-1)}}\) será el punto de partida para los cálculos de la capa \(\ell-1\).
Este proceso continúa hasta calcular los gradientes respecto a todos los parámetros \(W^{(\ell)}, b^{(\ell)}\) de la red (y respecto a la entrada \(x\) si fuera necesario).
Eficiencia: Backpropagation (modo reverso) es muy eficiente cuando tenemos muchas variables de entrada (parámetros) y una única salida escalar (la pérdida), ya que calcula todos los gradientes \(\frac{\partial L}{\partial \theta_i}\) con un coste computacional bajo (aproximadamente entre 2 y 3 veces el coste del forward pass). Por eso es el método estándar para entrenar redes neuronales.
(Nota: Librerías modernas como PyTorch, TensorFlow y JAX implementan AD de forma automática. El usuario define el paso hacia adelante (la estructura del modelo y el cálculo), y la librería se encarga de calcular los gradientes eficientemente usando backpropagation.)
Inicialización de parámetros
La forma en que inicializamos los pesos (\(W\)) y sesgos (\(b\)) de una red neuronal es sorprendentemente importante, sobre todo en redes profundas (con muchas capas). Una mala inicialización puede impedir que la red aprenda.
No podemos inicializar todos los pesos y sesgos a cero, ya que permanecerán iguales a lo largo del entrenamiento. En cambio inicializamos todos los sesgos a cero y los pesos de forma aleatoria. Pero esto tiene una sutileza.
El problema de la varianza grande o pequeña
La clave está en controlar la varianza de las salidas de cada capa (activaciones) y de los gradientes que fluyen hacia atrás. Aquí, \(D_\ell\) representa el ‘fan-in’ de la capa \(\ell\), es decir, el número de neuronas en la capa anterior que alimentan a una neurona de la capa \(\ell\).
Los pesos se inicializan de forma aleatoria, lo que significa que los pesos tendrán una varianza \(\sigma^2\), que podemos escoger y puede ser distinta para cada capa. Resulta que el valor de esta varianza es importante para el entrenamiento de redes profundas ya que:
Varianza demasiado grande (e.g., \(\sigma^2 = 5 / D_\ell\)): Las activaciones tienden a crecer exponencialmente capa tras capa (activaciones explosivas). Como la activación de cada capa se multiplica por la de la capa anterior, esto puede llevar a la saturación de las funciones de activación (como sigmoid o tanh, o incluso valores enormes con ReLU) y a gradientes gigantescos (gradientes explosivos), desestabilizando el entrenamiento.
Varianza demasiado pequeña (e.g., \(\sigma^2 = 0.1 / D_\ell\)): Las activaciones tienden a atenuarse exponencialmente capa tras capa (activaciones desvanecientes). Como la activación de cada capa se multiplica por la de la capa anterior, esto provoca que las señales (y sus gradientes) que llegan a las primeras capas sean minúsculas (gradientes desvanecientes), impidiendo que esas capas aprendan (como el gradiente es pequeño, los pesos no se actualizan).
La @fig_varianza_activaciones ilustra cómo evoluciona la varianza de las activaciones a través de las capas para diferentes inicializaciones en una red con activación ReLU:
Pequeña (\(\sigma^2 = 0.1 / D_\ell\)): La varianza de las activaciones se atenúa rápidamente hacia cero.
Xavier/Glorot (\(\sigma^2 = 1 / D_\ell\)): Mantiene la varianza más estable, pero está diseñado teóricamente para activaciones lineales o simétricas (como tanh). Con ReLU, tiende a reducir la varianza. (La versión original de Xavier promedia fan-in y fan-out: \(\sigma^2 = 2 / (D_\ell + D_{\ell+1})\)).
He (\(\sigma^2 = 2 / D_\ell\)): Diseñada específicamente para ReLU. Mantiene la varianza de las activaciones aproximadamente constante a lo largo de las capas. El factor 2 extra (comparado con Xavier) compensa el hecho de que ReLU anula ~la mitad de las entradas (las negativas), lo que de otro modo reduciría la varianza.
Grande (\(\sigma^2 = 5 / D_\ell\)): La varianza explota brutalmente, saliéndose de escala y haciendo inviable el aprendizaje.
Código
import numpy as npimport matplotlib.pyplot as pltnp.random.seed(42)# Parametersn_layers =20fan_in =100# D_\elln_samples =1000# Variance schemes to compareinit_schemes = {r"Pequeña ($0.1/D_\ell$)": 0.1,r"Xavier/Glorot ($1/D_\ell$)": 1.0,r"He ($2/D_\ell$)": 2.0,r"Grande ($5/D_\ell$)": 5.0,}colors = ["gray", "blue", "green", "red"]# Store resultsall_variances = {}# Simulate for each schemefor (label, scale), color inzip(init_schemes.items(), colors): variances = []# Input: standard normal x = np.random.randn(fan_in, n_samples) var = np.var(x) variances.append(var)for layer inrange(n_layers):# Initialize weights for this layer W = np.random.randn(fan_in, fan_in) * np.sqrt(scale / fan_in)# Linear + ReLU x = W @ x x = np.maximum(0, x) var = np.var(x) variances.append(var) all_variances[label] = variances# Plotplt.figure(figsize=(8, 5))for label, color inzip(init_schemes.keys(), colors): plt.plot(all_variances[label], label=label, color=color, linewidth=2)plt.xlabel("Capa (layer)")plt.ylabel("Varianza de activaciones")plt.title("Evolución de la varianza de activaciones según inicialización (ReLU)")plt.legend()plt.yscale('log')plt.grid(True, alpha=0.3)plt.tight_layout()plt.show()
Varianza de las activaciones según inicialización
Esto se puede entender porque la varianza de \(z^{(\ell)}\) está dada por:
donde \(\sigma_\ell^2\) es la varianza de los pesos de la capa \(\ell\). Esta relación también nos da una pista de cómo resolver el problema.
Derivación de la relación entre varianzas
Queremos encontrar una relación entre la varianza de las preactivaciones de la capa \(\ell\) y la varianza de las activaciones de la capa \(\ell - 1\). Sea \(\sigma_{\ell}^2\) la varianza de los pesos de la capa \(\ell\) y \(\sigma_{z^{(\ell)}}^2\) la varianza de las preactivaciones de la capa \(\ell\).
Para empezar, calculamos el valor esperado de la preactivación \(z^{(\ell)}\) en la capa \(\ell\):
Ahora bien, como \(ReLU\) corta a cero las entradas negativas, la varianza de las activaciones \(a^{(\ell)}\) es la mitad de la varianza de \(z^{(\ell)}\) (ya que la integral de \(x^2\) desde 0 hasta \(\infty\) es la mitad de la integral desde \(-\infty\) hasta \(\infty\)). Entonces:
Inicialización para el forward pass: \(\sigma_\ell^2 = 2/D_\ell\)
Para asegurar que la varianza de las activaciones \(a^{(\ell)}\) se mantenga aproximadamente constante durante el forward pass (al pasar de la capa \(\ell\) a \(\ell+1\)) cuando se usa la activación ReLU, Kaiming He et al. propusieron inicializar los pesos \(W^{(\ell)}\) muestreando de una distribución con media 0 y varianza:
donde \(D_\ell\) es el fan-in de la capa \(\ell\). Si reemplazamos esto en la @eq_varianza_z obtenemos \(\sigma_{z^{(\ell)}}^2 = \sigma_{z^{(\ell - 1)}}^2\), tal que la varianza es estable.
Inicialización para el backward pass: \(\sigma^2 = 2/D_{\ell+1}\)
De forma análoga, si analizamos el backward pass, para que la varianza de los gradientes \(\frac{\partial \mathcal{J}}{\partial z^{(\ell)}}\) se mantenga estable al retropropagar desde la capa \(\ell+1\) a la capa \(\ell\) (asumiendo activación ReLU), la condición requiere una varianza basada en el ‘fan-out’:
Notamos que la condición ideal para el forward pass (usa \(D_h\)) y para el backward pass (usa \(D_{h+1}\)) solo coinciden si \(D_h = D_{h+1}\). En la práctica, se puede usar la inicialización de He cuando la activación es ReLU, o un promedio \(\sigma_\ell^2 = \frac{4}{D_h + D_{h+1}}\).
Nota: Cuando la activación es tanh o sigmoide, la derivación de la relación entre las varianzas sugiere usar \(\sigma_\ell^2 = \frac{1}{D_h}\), que se llama inicialización de Xavier/Glorot.
Ejemplo de código de entrenamiento
Existen librerías que implementan todos los pasos de este entrenamiento. La más popular es tal vez PyTorch. Este curso no está enfocado en implementar el código ya que este cambia constantemente, y el simple uso de librerías se puede lograr leyendo la documentación, pidiendo ayuda en foros como StackOverflow o preguntándole a una IA. Sin embargo, damos algunos ejemplos, empezando por un ejemplo de cómo se implementaría lo visto hasta ahora usando PyTorch.
import torch, torch.nn as nnfrom torch.utils.data import TensorDataset, DataLoaderfrom torch.optim.lr_scheduler import StepLR# Definir el tamaño de la entrada, el tamaño de las capas ocultas y el tamaño de la salidaD_i, D_k, D_o =10, 40, 5# crear el modelo con dos capas ocultasmodel = nn.Sequential( nn.Linear(D_i, D_k), nn.ReLU(), nn.Linear(D_k, D_k), nn.ReLU(), nn.Linear(D_k, D_o))# Inicialización de Hedef ini_pesos(capa_in):ifisinstance(capa_in, nn.Linear): nn.init.kaiming_normal_(capa_in.weight, nonlinearity='relu') capa_in.bias.data.fill_(0.0)model.apply(ini_pesos)# Función de pérdida ECMcriterion = nn.MSELoss()# Construir el optimizador por SDG e inicializar la tasa de aprendizaje y el momentumoptimizer = torch.optim.SGD(model.parameters(), lr=0.1, momentum=0.9)# Objeto que reduce la tasa de aprendizaje por la mitad cada 10 épocasscheduler = StepLR(optimizer, step_size=10, gamma=0.5)# Crear 100 muestras aleatoriasx = torch.randn(100, D_i)y = torch.randn(100, D_o)dataloader = DataLoader(TensorDataset(x, y), batch_size=10, shuffle=True)# Entrenar el modelo, pasando 100 veces por los datosfor epoch inrange(100): epoch_loss =0.0for i, data inenumerate(dataloader): x_batch, y_batch = data# Poner a cero los gradientes optimizer.zero_grad()# Forward pass y_pred = model(x_batch)# Calcular la pérdida loss = criterion(y_pred, y_batch)# Backward pass loss.backward()# Actualizar los pesos optimizer.step()# Actualizar la pérdida epoch_loss += loss.item()# Imprimir la pérdidaprint(f"Epoch {epoch:5d}, Loss: {epoch_loss:.3f}")# Actualizar la tasa de aprendizaje scheduler.step()
Ejemplo aplicado de una red neuronal
En el libro hay un ejemplo para un problema de clasificación con 40 características. Recomiendo leerlo.
Otro ejemplo típico es la clasificación de dígitos escritos, el llamado MNIST. Consiste en un conjunto de imágenes en blanco y negro que representan dígitos escritos a mano. Se entrena una red para identificar a cuál dígito corresponde cada imagen.
Primero cargamos los datos. Note que necesitamos normalizarlos. En general las redes neuronales funcionan mejor si los datos tienen media cero y varianza unitaria. Por eso dividimos la media y la desviación estándar del conjunto de entrenamiento y usamos esas para normalizar el conjunto de entrenamiento y el de prueba.
Para lograrlo, primero cargamos los datos sin normalizar, les calculamos su media y varianza, y finalmente usamos transforms.Normalize para normalizar los datos.
import torchimport torch.nn as nnimport torch.optim as optimfrom torchvision import datasets, transformsimport matplotlib.pyplot as plt# Definir transformaciones para los datostransform = transforms.Compose([ transforms.ToTensor()])# Cargar datos de MNISTtrain_dataset = datasets.MNIST('./data', train=True, download=True, transform=transform)test_dataset = datasets.MNIST('./data', train=False, transform=transform)# Calcular media y desviación estándar del conjunto de entrenamientotrain_loader_temp = torch.utils.data.DataLoader(train_dataset, batch_size=len(train_dataset))data =next(iter(train_loader_temp))[0]mean = data.mean().item()std = data.std().item()# Aplicar normalización con los valores calculadostransform_norm = transforms.Compose([ transforms.ToTensor(), transforms.Normalize((mean,), (std,))])# Recargar los datasets con la normalización correctatrain_dataset = datasets.MNIST('./data', train=True, download=True, transform=transform_norm)test_dataset = datasets.MNIST('./data', train=False, transform=transform_norm)# Definir dataloaderstrain_loader = torch.utils.data.DataLoader(train_dataset, batch_size=64, shuffle=True)test_loader = torch.utils.data.DataLoader(test_dataset, batch_size=1000, shuffle=False)
Podemos graficar algunas de las imágenes del conjunto de entrenamiento.
# Graficar algunas imágenes del conjunto de entrenamientoimport matplotlib.pyplot as plt# Tomar 16 imágenes del conjunto de entrenamientoimages, labels =next(iter(train_loader))images = images[:16]labels = labels[:16]# Graficar las imágenesplt.figure(figsize=(10, 10))for i inrange(16): plt.subplot(4, 4, i +1) plt.imshow(images[i].numpy().squeeze(), cmap='gray') plt.title(labels[i].item()) plt.axis('off')plt.show()
Ahora vamos a entrenar un perceptrón multicapa para clasificar las imágenes.
# Definir una red neuronal MLP simpleclass MLP(nn.Module):def__init__(self):super(MLP, self).__init__()self.flatten = nn.Flatten()self.layers = nn.Sequential( nn.Linear(28*28, 128), nn.ReLU(), nn.Linear(128, 64), nn.ReLU(), nn.Linear(64, 10) )def forward(self, x): x =self.flatten(x)returnself.layers(x)# Inicializar modelo, función de pérdida y optimizadordevice = torch.device("cuda"if torch.cuda.is_available() else"cpu")# # Diagnóstico para disponibilidad de la GPU# print("Calculando en: ", device)# print("Dispositivo disponible: ", torch.cuda.is_available())# if torch.cuda.is_available():# print("Dispositivo actual: ", torch.cuda.current_device())# print("Nombre del dispositivo: ", torch.cuda.get_device_name(torch.cuda.current_device()))model = MLP().to(device)criterion = nn.CrossEntropyLoss()optimizer = optim.SGD(model.parameters(), lr=0.01, momentum=0.9)# Función para evaluar el modelodef evaluate(model, data_loader, criterion, device): model.eval() total_loss =0 correct =0 total =0with torch.no_grad():for data, target in data_loader: data, target = data.to(device), target.to(device) output = model(data) total_loss += criterion(output, target).item() * data.size(0) _, predicted = torch.max(output.data, 1) total += target.size(0) correct += (predicted == target).sum().item()return total_loss / total, correct / total# Entrenamientoepochs =10train_losses = []test_losses = []train_accuracies = []test_accuracies = []for epoch inrange(epochs): model.train() running_loss =0.0 correct =0 total =0for data, target in train_loader: data, target = data.to(device), target.to(device) optimizer.zero_grad() output = model(data) loss = criterion(output, target) loss.backward() optimizer.step() running_loss += loss.item() * data.size(0) _, predicted = torch.max(output.data, 1) total += target.size(0) correct += (predicted == target).sum().item() train_loss = running_loss /len(train_loader.dataset) train_accuracy = correct / total# Evaluar en conjunto de prueba test_loss, test_accuracy = evaluate(model, test_loader, criterion, device)# Guardar métricas train_losses.append(train_loss) test_losses.append(test_loss) train_accuracies.append(train_accuracy) test_accuracies.append(test_accuracy)# # Para ver el progreso del entrenamiento# # comentarlo a la hora de producir el sitio web# print(f'Epoch {epoch+1}/{epochs}: '# f'Train Loss: {train_loss:.4f}, Train Acc: {train_accuracy:.4f}, '# f'Test Loss: {test_loss:.4f}, Test Acc: {test_accuracy:.4f}')
Después de entrenar, graficamos las pérdidas y la precisión en el conjunto de entrenamiento y prueba.
# Graficar resultadosplt.figure(figsize=(10, 6))plt.subplot(1, 2, 1)plt.plot(train_losses, label='Pérdida en entrenamiento')plt.plot(test_losses, label='Pérdida en prueba')plt.xlabel('Época')plt.ylabel('Pérdida')plt.legend()plt.title('Pérdida vs. Época')plt.subplot(1, 2, 2)plt.plot(train_accuracies, label='Precisión en entrenamiento')plt.plot(test_accuracies, label='Precisión en prueba')plt.xlabel('Época')plt.ylabel('Precisión')plt.legend()plt.title('Precisión vs. Época')plt.tight_layout()plt.show()